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Abstract. We present an improved prediction of the nonlinear perturbation theory (PT) 
via the Lagrangian picture, which was originally proposed by Matsubara (2008). Based on 
the relations between the power spectrum in standard PT and that in Lagrangian PT, we 
derive analytic expressions for the power spectrum in Lagrangian PT up to 2-loop order in 
both real and redshift spaces. Comparing the improved prediction of Lagrangian PT with 
A^-body simulations in real space, we find that the 2-loop corrections can extend the valid 
range of wave numbers where we can predict the power spectrum within 1 % accuracy by a 
factor of 1.0 {z = 0.5), 1.3 (1), 1.6 (2) and 1.8 (3) vied with 1-loop Lagrangian PT results. 
On the other hand, in all redshift ranges, the higher-order corrections are shown to be less 
significant on the two-point correlation functions around the baryon acoustic peak, because 
the 1-loop Lagrangian PT is already accurate enough to explain the nonlinearity on those 
scales in A^-body simulations. 
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1 Introduction 

In the forthcoming era of precision cosmology, accurate and precise predictions of cosmological 
theories play crucial roles not only to estimate the cosmological parameters but also to detect 
a tiny signal from the observational data. The large-scale structure of the universe is a 
powerful probe of cosmology, and recently the baryon acoustic oscillations (BAOs) imprinted 
on the large-scale structure [1-3] have attracted much attention to constrain the nature of 
dark energy. 

There are many ongoing and planned redshift surveys aiming at clarifying the nature 
of the dark energy through a precision measurement of BAOs, including the Sloan Digi- 
tal Sky Baryon Oscillation Spectroscopic Survey (BOSS) [8], Hobby-Eberly Telescope Dark 
Energy Experiment (HETDEX) [9], BigBOSS [10], Subaru Measurement of Images and Red- 
shifts (SuMIRe) [12], Wide-Field Infrared Survey Telescope (WFIRST) [11], Euchd [13]. To 
measure the characteristic scales of BAO precisely enough for investigating the dark energy, 
accurate theoretical template for the power spectrum and correlation function on BAO scales 
is quite essential. 

Although the nonlinear gravitational clustering of galaxies on BAO scales is rather 
moderate, the linear theory of density perturbations [14] is inadequate to quantitatively 
describe the BAOs [15, 16]. The standard perturbation theory (SPT) provides a systematic 
way to investigate the gravitational clustering of dark matter in nonlinear regime [17-23, 25, 
26]. Even though the SPT describes the nonlinearity fairly well at sufficiently high redshifts 
[27], it is still insufficient to describe the clustering on BAO scales at observationally relevant 
redshifts, z = 0-3. To overcome such situation, various methods to improve the SPT are 
proposed, partially resumming infinite series of higher-order perturbations. The renormalized 
perturbation theory (RPT) [28-30] is a representative pioneering work, and since it appeared, 
various approaches have been subsequently proposed [31-38]. The quantitative predictions of 
those approaches are compared with A^-body simulations, and actually shown to excellently 
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reproduce the A^-body results quite well in a certain range of power spectrum and two-point 
correlation function. The accuracy and precision of the predictions actually depend on the 
treatment of each method, and a comprehensive test of various methods is examined in 
ref. [46]. 

Most of the above newly developed perturbation theories predict the nonlinear gravi- 
tational clustering of dark matter in real space. However, the observable quantity in actual 
redshift surveys is the clustering of galaxies in redshift space. There are two ingredients to 
be added to the perturbation theory in order to compare the theoretical predictions with 
observations: redshift-space distortions [47-49] and biasing [51, 52]. While the redshift dis- 
tortions induce the clustering anisotropics, the galaxy biasing can produce a scale-dependent 
enhancement (or suppression) in the clustering amplitude. Hence, coupled with the effect of 
gravitational clustering, both the amplitude and shape of the power spectrum and two-point 
correlation function would be significantly changed, and an accurate modeling of BAOs that 
ensures a percent-level precision seems non-trivial. 

The resummed PT via the Lagrangian picture provides a unique framework to include 
both the ingredients into the fundamental formulation, and does not have any phenomeno- 
logical parameter of dynamics once the biasing model is fixed [37, 53, 54] (but see [50]). 
In contrast to the SPT formulated as Eulerian PT, the resummed PT by refs. [37, 53, 54] 
is based on the framework of the Lagrangian perturbation theory (LPT) [55-63], in which 
the displacements of fluid element are treated as dynamical variables. It is known that the 
first-order LPT reproduces the classic Zel'dovich approximation [64]. In both Eulerian and 
Lagrangian PTs, the small variables in the perturbative expansion correspond to the initial 
density fluctuations. Therefore, when the theoretical predictions of LPT are made in the 
Eulerian space with the full expansion of perturbations, LPT should give mathematically 
equivalent results to the SPT expansion [54]. In the Lagrangian resummation method, corre- 
lations of the displacement field at zero-lag are kept in the exponent, and other contributions 
to the power spectrum are expanded as usual in Eulerian space. This partial expansion cor- 
responds to the resummation of a class of infinite higher-order perturbations. As a result, 
an exponential prefactor naturally arises in the expression of the nonlinear power spectrum 
[37]. The exponential prefactor plays an important role in describing the nonlinear smearing 
effects of BAO in both real and redshift spaces. 

In previous works, the Lagrangian resummation method has been investigated only up 
to 1-loop order in perturbations. The purpose of this paper is to generalize the previous work 
and to carry out 2-loop calculations of the same resummation method. We derive analytic 
expressions of the 2-loop results in both real and redshift spaces. To see how well the 2-loop 
corrections improve the previous 1-loop results, we focus on the power spectrum of dark 
matter in real space, and compare the predictions with numerical simulations. Although 
we leave a quantitative analysis in redshift space for future work, the improvement of the 
predictions shown in this paper is an important step, and would lead to a strong motivation 
to actually evaluate the effects of redshift-space distortions and the biasing via the LPT up 
to 2-loop order. 

This paper is organized as follows. In section 2, useful equations and concepts, which 
are mostly derived in previous work, are summarized. In section 3, our main results of the 
2-loop calculations both in real and redshift spaces are presented. In section 4, the analytic 
expressions are compared with A^-body simulations of dark matter. In section 5, our results 
and conclusions are summarized. 
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2 Lagrangian Perturbation Theory and Standard Perturbation Theory 



In this section, we summarize some concepts and equations of Lagrangian resummations 
method, which we use in this paper. Readers find derivations of the following equations in 
refs. [37, 53]. 

According to the result of 1-loop resummation via LPT, the following identity is con- 
firmed by straightforward calculations of 1-loop SPT and LPT: 



-Plpt(^) = exp 



(2.1) 



where PhPTik) is the nonlinear power spectrum predicted by the 1-loop resummation via 
LPT, Phik) is the linear power spectrum, and P^p!^^{k) is the 1-loop contributions (without 
linear contribution) to the power spectrum predicted by the SPT. Note that the time depen- 
dence of the power spectrum Plpt is implicitly encoded in the linearly extrapolated power 
spectrum Phik) evaluated at the redshift of our interest. 

As we mentioned in section 1, the displacement field is the fundamental variable in 
LPT, and is treated as a small perturbative quantity. In computing the power spectrum, we 
need to evaluate the correlation of displacement fields at different positions in the exponent. 
The resummation method proposed by refs. [37, 53] gives a specific but efficient recipe to 
treat this. That is, the correlation of the displacement fields at a single position is kept in 
the exponent, while the correlation among components of separated displacement vectors, 
which is expected to be small compared to zero-lag correlation as long as we consider a large 
separation, is expanded from the exponent. The exponential prefactor in eq. (2.1) indeed 
represents the zero-lag correlation up to the 1-loop level. It is easily checked that expanding 
the exponential prefactor and truncating the series up to the 1-loop order exactly recovers 
the power spectrum of 1-loop SPT. 

Note that in the language of RPT [28, 29], the resummation treatment by refs. [37, 53] 
corresponds to a partial vertex renormalization of the SPT expansion [37], and the expo- 
nential prefactor in eq. (2.1) contains the information beyond the 1-loop SPT. Further, the 
resultant functional form of the exponential prefactor quite resembles the renormalized prop- 
agator at high-A: limit [29], and is shown to play a physically important role in describing 
the nonlinear smearing effect of BAOs. Nevertheless, because of a rather strong damping, 
the resultant predictions for the power spectrum ceases to work well at higher-A; modes, 
where the non- linear enhancement of the power spectrum amplitude becomes significant. 
Thus, similar to RPT and other improved PT treatment, the accuracy of the Lagrangian 
resummation treatment will be limited to a specific range k < fedampi where fcdamp is the 
characteristic scale of the exponential damping factor. But, this characteristic scale depends 
on the order of loop corrections included in the power spectrum expression, and including 
the next-to-leading order will definitely improve the predictions applicable to some higher-A: 
modes. 

The main purpose of this paper is to extend the 1-loop expression (2.1) to a general 
relation between the power spectrum in SPT and that in LPT for arbitrary loop-order, and 
is to explicitly derive the 2-loop expression. To do so, we first notice that the exponen- 
tial prefactor corresponds to the square of a Fourier transform of the one-point probability 
distribution function of the displacement field, (e"*'" '^)^, where ^{q) = x{q) — q is the 
displacement field, x and q are respectively Eulerian and Lagrangian coordinates of a fluid 
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element [53, 54]. According to the cumulant expansion theorem [65], this factor reduces to 



exp 



n=l 



(2n)! 



(2.2) 



where 



A) 



(2n) 



■ 2n 

n 

.1=1 



(27r)3 



2n 



(2.3) 



and Ci^...i^^ is a 2?i-order cumulant of the displacement field in Fourier space defined by 



'^*i(Pi)---'^*2.(P: 



'2n, 



The square of eq. (2.2) is naturally factorized out in LPT for the nonlinear power 
spectrum. On the other hand, the LPT and SPT should give the same result when the 
perturbations are completely expanded in Eulerian space, since the fundamental equations 
of motion, and small perturbations of the expansion are the same in both perturbation 
schemes [54]. Thus, when we factorize out the square of eq. (2.2) from SPT series of the 
power spectrum, we have the following identity: 



P{k) = exp 



■2E 

n=l 



(2n) 

n---«2»i 



(2n)! 

OO 



m=l 



1 + 



OO 

1=1 



n=l 



(2n)! 



ki2„ A{2n) 



(2.5) 



where P^}^°^{k) is an m-loop contribution to the power spectrum in SPT. The indices 
J) ^1) ^2 etc. correspond to spatial components and run over 1, 2, 3. Repeated appearance of 
zi,i2, . . . implicitly means those indices should be summed over. 

As we mentioned, in the Lagrangian resummation method, the exponential prefactor is 
kept unexpanded, and perturbative truncations are made in the exponent. The quantity in 
the last bracket is normally expanded and truncated as usual. In the A^-loop approximation 
of the Lagrangian resummation, the series in the exponent is expanded up to 0{Pi,)^ , and 
the remaining factor is expanded up to ©(Pl)^"*"^- The case of = 1 just reduces to 
eq. (2.1). 

Li this paper, we generalize eq. (2.1) to the case of 2-loop approximation, N = 2. From 
the consideration above, we do not need to calculate LPT from the beginning. Instead, we 
first calculate P^}^°'^{k) for m = 1,2, and use eq. (2.5) truncated at the N = 2 order to 
obtain 2-loop predictions of the Lagrangian resummation method. 

We need the polyspectra Cji...i2„ of only n = 1, which are given by 



Cijip, -p) = Cij{p) 



4^^(P) + C^f(p) + Cif)(p) + C^f)(p), 



(2.6) 
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where 



cfp' 



C\f\p) = \j ^L't\p',p-p')Lf\p',p-p')PUp')PL{\p-p'\), 



(2.7) 
(2.8) 



4=^)(p) = cjf)(p) 



h\'\p)PUp) 



^Lf\p,-p',p')PUp') 



(2.9) 



up to ©(Pl)^- The perturbative kernels in LPT up to third order are given by [61, 62] 



Lf\pi,P2) 



Lf""^ (Pl,P2,P3) 









3 ki 


1 - 


fPl-P2 




V PlP2 


5 ki 


1 - 


fPl-P2 


7fc2 


V PlP2 



(2.10) 
(2.11) 



(Pl + P2) ■ P 3 

\Pl +P2\P3 



1 2 



3A;2 



1 



Pl-P2\" ^ ^ (Pl •P 2)(P2 •P3)(P3 -Pi) 



P1P2 

+ fc X T(pi,p2,P3), 

-(n) 



(2.12) 



where fc = pi + • • • +Pn for L ■ , and a vector T represents a transverse part whose expression 
is not needed in the foUowing apphcation. The third-order kernel L^^^ in eq. (2.9) is a 
symmetrized one in terms of their arguments: 



4^^(Pl,P2,P3) 



-^^f''^(Pi,P2,P3) + perm. 



(2.13) 



Next consider the power spectrum in redshift space. In the Lagrangian picture, the 
displacement field in redshift space, is distorted by the peculiar velocities, and mapped 
from that in real space as 



= * + 



z * 

IT' 



(2.14) 



where H is is the Hubble parameter, the unit vector z indicates the line-of-sight direction 
and ^ = d^/dt. Using the fact that the time-dependence of the perturbative kernel in real 
space is approximately described by oc with D being the linear growth factor, the 
displacement field of each perturbation order in redshift space is given by 



(2.15) 



where / = din D/ din a is the logarithmic derivative of the linear growth factor. Then, the 
mappings from real space to redshift space of the order-by-order cumulants are given by 



c: 



s{nm) _ p(Ti) p(m)^(mn) 
ij ~ ^ik ^jl ^kl ' 



(2.16) 
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where 



R 



(n) 



6ij + nfziZj. 



(2.17) 



Thanks to the hnear mapping (2.15), the perturbation calculation of the displacement field 
is rather straightforward in redshift space, provided the perturbative solutions in real-space. 
Nevertheless, the resultant form of the redshift-space power spectrum is rather complicated 
because of the coupling between velocity and density fields. The general expressions for 
redshift-space power spectrum is formally written as 



Ps{k) = exp 



CXJ J 

n=l 



(2n)! 



m-loop 
sSPT 



m=l 

oo 



\ ^ (2n)! 

1=1 \ n=l ^ ' 



(k) 



, (2.18) 



where -Ps™pt°^(^) m-loop contribution to the redshift-space power spectrum in SPT. 



The quantity ^i|^."2„ 



is defined by 



A' 



(2n) 



■ 2n 

n 



(2vr)= 



6l 



2n 



(2.19) 



where Cl^...j^^^ is defined by a similar mapping as that of eq. (2.16). 

Although the above result involves many additional terms compared to the real-space 
power spectrum, several important terms in characterizing redshift distortions is identified in 
eq. (2.18). One is the enhancement factor, (l + //i^)^, known as the linear Kaiser factor [47]. 
The corrections to the exponential prefactor give an additional damping of the power spec- 
trum amplitude, which roughly corresponds to the Finger-of-God effect driven by the random 
motion of fiuid element. In section 3.2, we will explicitly see the corrections characterizing 
the redshift distortion up to the 2-loop order. 



3 Analytic expressions of the 2-Ioop power spectrum in Lagrangian re- 
summation method 

Based on the relation between the power spectrum in the SPT and that in the LPT, the 
2-loop LPT power spectrum can be expressed as a combination of 1- and 2-loop SPT power 
spectra. In this section, we derive analytic expressions of the real- and redshift-space power 
spectra up to the 2-loop order in LPT. 

3.1 Real space 

Let us first evaluate the power spectrum in real space. The only task is to calculate the terms 
in the exponent of eq. (2.5) up to the demanding order in Ph{k). Up to the second order in 
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Phik), the exponent of eq. (2.5) is expanded as ^ 

L-A-. 

J ^3 12 



- (2n)! 



6^ 



2 r 



^(11) ^ ^(22) ^ 2^(13) 



The quantity A^"^^^ is calculated from the polyspectra C^-"^'' as 



(3.1) 



(3.2) 



where we define C|" (p) = C°'^{p)piPj and p.j = Pi/p. The explicit expressions for C"'^ are 
presented in eqs. (A9)-(A11) of ref. [37]. In the above, we used the fact that f d^pj {'27r)^piPjf[p) 
(5ij/(67r^) J p'^dpf{p) for an arbitrary function f{p). 

Based on the perturbative kernels in LPT defined by eqs. (2.10)-(2.13), the 1-loop 
correction term, A^^^\ is explicitly written as [37] 



dp Pl{p)- 



(3.3) 



The next-to-leading order contributions, A^'^'^^ and A^^'^\ are also given by 



^(22)= / p^dpC^^^\p) 



p dp- 



1 f d^pi 



2 J (27r) 



Lf\pi,p-pi) LX>{pi,p-pi) PUpi)Pl{\p-Pi\) 



(2) 



^(13) 



3927r2 Jo 

p^dp C^^^^p) 



PiP2dpidp2Ph{pi)Pi.{P2)K \{p\Ip2 +P2/Pl)/'A 



(3.4) 



p'dp\5,,Lf\p)Pi^{p) I ^Lf\p,-p^,p^)PUpi) 

PlP2dpidp2PL{Pl)PL{p2)K [(Pl/P2 +P2M)/2] , 



847r2 



where the function K{y) is defined by 



y-1 



y+l 



(3.5) 



(3.6) 



We note that the function K{y) is always positive in a range y > 0. Also, the 2-loop 
corrections, 

^(13) and ^(22) 

are shown to be positive. 



-^ijim 'ioes not have a C'(Pl)^ connected cumulant. 
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Collecting the 2-loop corrections given above, the power spectrum up to 2-loop order in 
LPT is expressed as 



2-loop 



LPT 



(k) 



exp 



^(11)^^(22)^2^(13)^ 



PL(A:)+p^p'rw+^spr(^) 



SPT 
„2 



+ PUk) I ^ [a^''^ + A^''^ + 2A^''^) + 



k^ 



(3.7) 



where P^^?^^{k) and fgpx°''(^) s-i'g 1" ^iid 2-loop contributions to the real space power 
spectrum in SPT respectively. 

3.2 Redshift space 

Next consider the power spectrum in observable redshift space. In similar manner to the real 
space, we start with eq. (2.18) by expanding the exponent up to the second order in Pi^{k): 



-2E 

n=l 



' ' ' ^i2n /is(2n) 

{2n)\ 



s(2) 



-kikjA-j 



^2 ^i^j^l^mA, 



s(4) 
ijlm 



+ 



hh r I 4«(22) , /is(13) , .s(31) 



(3. 



The quantity A^^°^^^ represents a momentum integration of a polyspectrum in redshift space. 
Based on the relation between the polyspectra in real and redshift spaces, we obtain 



(27r)3 *J 

{6ij + {a + (3)fz,zj + af3fziz,} A^^^^l 



Taking an inner product of A^^"^^^ with kikj , we have 



k,k,Af^^ =k'V^'^^\f,fi)A^''''\ 



(3.9) 



(3.10) 



The quantity fi = z ■ k/k is the directional cosine of the wave vector k with respect to the 
line of sight. The function 'D^'^^\i^i, f) which describes the redshift distortion is defined by 



V^^^\^^, f) = p(^")(/i, /) ^ 1 + (a + /3 + aPf)ffi\ 



(3.11) 



From these relations, we can get the non-linear power spectrum in redshift space up to 2-loop 
order. The result is 

+ p(22)^(22) ^ 2P(13)^(13)| 



-^sLPt'^C^) = exp 



32-loop / 



+ (l + //iyPL(fe) 
k^ 



67r2 



_^ p(22)^(22) _^ 2P(^2)^(^^) 



fc2 



6^2 



(3.12) 
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where -Psspx''(^) and -Ps^px''(fc) are 1- and 2-loop contributions to the redshift-space power 
spectrum in SPT, respectively. 

In eq. (3.12), the terms proportional to / /i^ represent the corrections characterizing 
the redshift distortions. These corrections appear not only in the spatial-correlation terms, 
which have been expanded from the exponent, but also in the exponential prefactor. Thus, 
the resultant redshift-space power spectrum naturally possesses an additional damping effect 
responsible for the Finger-of-God (FoG) effect. Note that a naive perturbative treatment 
based on the SPT usually fails to describe the FoG effect, which should be phenomenologically 
introduced by hand (but see [50]). By contrast, in LPT, the redshift distortion is described 
by a simple linear mapping from real space, (2.14) or (2.15), and with the resummation 
method, the power spectrum damping by the FoG effect is naturally incorporated into the 
redshift-space power spectrum in a non-perturbative manner. 

While the resultant expression for redshift-space power spectrum seems to contain phys- 
ically relevant effects of non-linear redshift distortions, an explicit evaluation of eq. (3.12) 
needs a bit technical aspect. To be specific, a difficult part is the term -Pss'pT^(^)- involves 
six-dimensional integral, and is characterized as a function of k and directional cosine fi. No 
one has so far tried to evaluate it, and at present, the explicit calculation of -Pslpt''(^) seems 
infeasible unless we develop an efficient technique to compute ^'s^pt''(^)- Hence, we leave 
the discussion on the numerical evaluation of the redshift-space power spectrum for future 
work, and we hereafter focus on the quantitative prediction of real-space power spectrum. 

3.3 Evaluations of the 1- and 2-loop power spectra in SPT 

The expression of real-space power spectrum given in eq. (3.7) involves the 1- and 2-loop 
power spectra in SPT, and we need to directly evaluate them for a quantitative prediction 
of the power spectrum. Here, we briefly explain how to compute these two contributions. 
The explicit expressions for the 1- and 2-loop corrections in SPT are decomposed into several 
pieces as 

-^SPT ^i^) ~ -^SPt(^) ~^ -PsPt(^)' 

-^SPT '^(^) ~ Pspt(^) + -^sptI^) + -fsPT(^)' (3.13) 

where the quantities -Pgp^^ imply the ensemble averages obtained from the m-th and n-th 
order perturbative solutions of SPT. They are expressed as a collection of multi-dimensional 
integrals involving the perturbative kernel of the higher-order solutions (e.g., [6, 24, 46]): 

^ / (03 {^'^'Hq,k-q)y PUq)PL{\k-q\), (3.14) 

6Pl(A;) I ■0^T^'Hk,q,-q)PLiq), (3.15) 
^^^^^^ {/ (03-^^'^(^'P'-p)^l(p)} 

+ 6 J {T^'\p,q,k-p-q)y PUp)PLiq)PLi\k-p-q\), (3.16) 

2^1 ^j0^T^^\p,k-p)T^''\p,q,-q,k-p)PUp)PLiq)PL{\k-p\), (3.17) 



-^spt(^) ~ 

Pgp^(/C) = 

-^spt(^) ~ 
-^spt(^) ~ 
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Table 1. Parameters of A^-body simulations. 



Name -^^box # of particles Zinit # of runs 



Main 


1,000/1" 




5123 


31 


30 


High (Lll-Nll) 


2, 048/1" 


"^Mpc 


2048^ 


99 


1 


(L12-N11) 


4, 096/1" 


"^Mpc 


2048^ 


99 


1 



Pg'^ik) = 30Pl(A;) I ^^J'^''\p,q,k,-p,-q) Pi.{p) Pi.{q) . (3.18) 

Here, J-^'^\ki, • • • , is the symmetrized kernel of the n-th order solutions for the density 
field (see [6, 26] for derivation and explicit expressions). 

Note that the expressions of the 1-loop power spectra -PgpT°^ can be further reduced 
to the one-dimensional and two-dimensional integral for P^^^) and respectively (e.g., 

[21-23]). For the explicit computation of power spectrum in next section, we use the method 
of Gaussian quadratures for numerical integration of 1-loop power spectra. On the other 
hand, for the 2-loop power spectra, except for the first term in P^^^\ the integrals are 
no longer simplified, and we directly evaluate the five-dimensional integration (due to the 
symmetry around the vector fc, one of the azimuthal angles is trivially integrated). In this 
paper, following the approach by [6], we adopt the Monte-Carlo technique of quasi-random 
sampling using the library Cuba^. The symmetrized kernels for perturbative solutions higher 
than the third order are too complicated to express analytically, and we numerically generated 
the symmetrized kernels using the recursion relation of perturbative solutions (e.g., [6, 20, 
26, 28]). 

4 Comparison with A^-body simulations in real space 

In this section, in order to quantify the effect of the next-to-leading order corrections in LPT, 
we explicitly compute the 2-loop corrections, and compare the analytic predictions with A^- 
body simulations. We adopt a cosmological model with the WMAP Syr [4] parameters 
= 0.279, l^A = 0.721, Jib = 0.046, h = 0.701, = 0.96 and cts = 0.817, and the linear 
power spectrum Ph{k) is calculated from the output of the CAMB code [5]. 

4.1 A^-body simulations 

To compare the LPT results with A^-body simulations, we use the subset of the A^-body data 
presented in [6]. The data were created by a public A^-body code GADGET2 [42] with cubic 
boxes of side length 1 /i~^Gpc, 512'^ particles. The initial conditions were generated by the 
2LPT code [43] at 

^init — 31, and the output data of the 30 independent realizations were 
stored at redshifts z = 3, 2, 1, and 0.5. In addition to the main simulations, we also use the 
results of high-resolution simulation data presented in ref . [44] , which were similarly created 
by GADGET2 with the initial condition generator, 2LPT. We summarize our sets of A^-body 
simulation in Table 1. 

The power spectra of these simulations were computed by combining several output 
results with different box sizes, and for the purpose to probe BAOs, we specifically use the 
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results of A^-body runs called Lll-Nll and L12-N11 in ref. [44], whose box sizes are 2, 048 and 
4, 096 Mpc, respectively. The simulations werG 6Volv6ci from the redsliift ^init — 

99 and 

the output results were stored at z = 3 and 1. Although the realization of these simulations 
is only one, each simulation contains 2,048'^ particles, and we preferentially use the Lll-Nll 
run to measure the small-scale power spectrum. 

We measure the matter power spectrum adopting the same treatment as described in 
ref. [40]. That is, we calculate the power spectrum from the Fourier transform of the density 
field assigned on the 1, 024^ grids with the Cloud- in-Cells interpolation. To reduce the effect 
of finite-mode sampling advocated by ref. [39] , the resultant power spectrum is multiplied by 
the ratio P\in{k) / P^^ [k) at k < 0.1/iMpc~^, where the quantity P^^{k) is calculated from 
the SPT up to the third-order in density field, and -Piin(^) is the input linear power spectrum 
extrapolated to a given output redshift. In computing P^'^{k), we use the Gaussian-sampled 
density field used to generate the initial condition of each A^-body run (see refs. [39, 40] 
in detail). As for the estimation of two-point correlation function, we use the grid-based 
calculation using the Fast Fourier Transformation (FFT) [6]. Similar to the power spectrum 
analysis, we first compute the square of the density field on each grid of Fourier space. 
Then, applying the inverse Fourier transformation, we take the average over realization, and 
obtain the two-point correlation function. Note that the finite-mode sampling also affects the 
calculation of the two-point correlation function. We correct it by subtracting and adding 
the extrapolated linear density field as, ^(r) — S,Vmi^) + S,hni^)i where S^nn is the correlation 
function estimated from the Gaussian density field, and is the linear theory prediction of 
two-point correlation function. 

4.2 Power spectrum 

Figure 1 shows the ratio of power spectrum to a smooth reference spectrum, P{k)/ Pnvj{k), 
where the function P^^{k) is the linear power spectrum calculated from the smooth transfer 
function neglecting the BAO feature in ref. [7]. From top to bottom, the results at redshifts 
z = 3, 2, 1 and 0.5 are shown. In each panel, the power spectra from linear theory (black 
dotted), 1-loop SPT (blue dashed), 2-loop SPT (blue sohd), 1-loop LPT (red dashed) and 
2-loop LPT (red solid) are plotted, and are compared with A^-body simulations. Results of 
the main and high-resolution simulations are respectively shown as black and gray symbols 
with error bars. Note that owing to 30 independent realizations and the correction of the 
finite- mode sampling by ref. [39], the scatter of the A^-body results is rather reduced, and 
the size of each error bar becomes hard to see visually. Overall, the predictions of the 2-loop 
LPT show better agreement with A^-body simulations. The range of the agreement is wider 
than that of the 1-loop LPT at all redshifts, indicating that the 2-loop corrections in LPT 
can give an important contribution to the nonlinear enhancement of the power spectrum at 
high k. On the other hand, the 2-loop correction in SPT slightly reduces the amplitude of 
power spectrum, and the predictions of 2-loop SPT tend to reproduce the A^-body results 
quite well compared to the 1-loop SPT. However, a closer look at the behaviors on small 
scales reveals that the predictions overestimate the A^-body results at high redshift, and then 
turn to slightly underestimate at low redshift. Note that the power spectra obtained from 
the main and high-resolution simulations show somewhat different behaviors at z = 3. At 
k > 0.2/iMpc-\ the result of main simulations gets smaller power, and the discrepancy 
between simulation and the 2-loop predictions apparently manifests. Presumably, this would 
be explained by a transient phenomenon due to the lack of number of particles. Low resolution 
simulations with a small number of particles usually suffer from the discreteness effect, and the 



- 11 - 






' ' ' ' 1 


/ ^ 




z=1.0 


1 F 








c 
D_ 






\ 






qT 




r \ ■ \ 




. . . . 'i 


. . : . , V . . . , A. . . ■ 



0.1 



0.2 
k [h/Mpc] 



0.3 



0.4 




0.4 



Figure 1. Ratio of power spectrum to smoothed reference spectrum, P{k)/ Pnw{k), given at redslrifts 
z ~ 3 (top left), 2 (top right), 1 (bottom left) and 0.5 (bottom right). The power spectra from 
linear theory (black dotted), 1-loop SPT (blue dashed), 2-loop SPT (blue sohd), 1-loop LPT (red 
dashed) and 2-loop LPT (red solid) are compared with A^-body simulations. Results of the main 
and high-resolution simulations are shown as black and gray symbols with error bars, respectively. 
The reference power spectrum Pnw(fc) is calculated from the no- wiggle formula of the linear transfer 
function in ref. [7]. 



A^-body data on the scales below the particle mean separation are not reliable particularly at 
an early time. Later, as the gravitational clustering develops, the non-linear mode transfer is 
expected to eventually dominate the initial power, and the small-scale structure tend to catch 
up the actual nonlinear growth (e.g., [45]). Although this qualitative picture is specifically 
relevant for the small-scales clustering, the apparent lack of initial power in low-resolution 
simulations may be influential even on relatively large scales, k > 0.2/iMpc~^, as shown in 
detail by ref. [44] (see figure 3 of their paper). It can affect the estimation of the validity 
range of PT results, and we need to take care of the systematics. 

In figure 2, to investigate the agreement in more quantitative ways, we show the dif- 
ferences of the power spectrum between simulation and PT prediction normalized by the 
smooth reference spectrum, i.e., [-P/v-body — -PpT]/-fnw(^)- Here, the differences between the 
main simulation and predictions are mainly shown, but we also plot the case with high- 
resolution simulation for 2-loop LPT ( main vs 2-loop SPT: blue, main vs 1-loop LPT: green, 
main vs 2-loop LPT: red, high-resolution vs 2-loop LPT: light red). In each panel of figure 
2, the green and red vertical arrows respectively represent the characteristic wave numbers 
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Figure 2. Difference between tlie results of main TV-body simulation and PT (2-loop SPT: blue, 1- 
loop LPT: green, 2-loop LPT: red) divided by the reference spectrum, [Pn -body (k) ~ Pprik)]/ Pnw{k)- 
We also show the difference between the result of high-resolution simulation and 2-loop LPT with 
light red symbols. In each panel, green and red vertical arrows represent the wave number k^^^°°^ /2 

and fc^j'°°^/2, which are respectively estimated from the exponential pre-factor in the 1- and 2-loop 
expressions for real-space power spectrum in LPT [see eqs.(4.1) and (4.2)]. 
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(4.1) 
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These wave numbers indicate the characteristic damping scales of the power spectrum in 
LPT, which are natural extension of the definition of the nonlinear scale, k^i, in ref. [37] [see 
eq. (38) of this paper]. 

Compared to the predictions of l-loop LPT, the 2-loop corrections in LPT extend the 
range of agreement with 1% precision by a factor of 1.0 {z = 0.5), 1.3 (1), 1.6 (2) and 1.8 (3). 
The ongoing (e.g., HETDEX [9]) and proposed surveys (e.g., BigBOSS [10], WFIRST [11], 
and SuMIRe [12]) plan to precisely measure the BAOs at z = 1-2, and thus the improvement 
with the 2-loop correction at high-z is practically important and encouraging for an accurate 
determination of acoustic scales. Figure 2 also implies that LPT has a better convergence 
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property for the higher-order corrections, since the range of agreement becomes monotonicahy 
wider as we include the higher-order corrections. Apparently, at lower redshifts, the range 
of agreement of the 2-loop SPT becomes comparable to that of the 2-loop LPT. As seen in 
figure 1, however, the 2-loop corrections in SPT can give a negative contribution at large k, 
and it is not guaranteed that the inclusion of higher-order corrections in SPT monotonically 
improves the prediction. In this respect, the agreement in SPT shown at lower redshifts may 
be regarded as accidental one. As the non-linearity significantly develops, even the 2-loop 
LRT is insufficient to improve the validity range of prediction. Although we did not store 
the simulation data at z = 0, it is expected that the validity range of 2-loop LPT at z = 
becomes rather comparable to that of 1-loop LPT, and the 1% precision would be achieved 
only at k< 0.7/i~^Mpc (see ref. [46]). 

Finally, it is worth mentioning a relation between the validity range of LPT and the 
exponential damping scale. Ref. [37] pointed out that the damping scale is considered as a 
criteria for the validity range of the 1-loop LPT, and we can see that the kl^^""^ /2 becomes 
a plausible indicator of the agreement with A^-body simulations. On the other hand, figure 
2 shows that such a criteria does not simply work for the 2-loop LPT. Though the damping 
scale of the 2-loop LPT is always smaller than that of the 1-loop LPT, due to the higher-order 
corrections in the exponential prefactor, the validity range of the 2-loop LPT becomes wider 
than the scale indicated by /c^J^°°^/2 or fc^j'°°^/2. In this sense, the damping scale k'^^°°^ /2 
may be regarded as the scale where the prediction of 2-loop LPT starts to deviate from that 
of the 2-loop SPT. 

4.3 Correlation function 

The two-point correlation function ^(r) is given by Fourier transform of the power spectrum: 



Figure 3 shows the two-point correlation functions around the baryon acoustic peak at red- 
shifts z = 3 (top left), 2 (top right), 1 (bottom left), and 0.5 (bottom right). In each panel, 
the predictions calculated from linear theory (black dotted with error bars), 1-loop LPT 
(green dashed) and 2-loop LPT (red solid) are plotted. Also, figure 4 shows the fractional 
difference between the A^-body and Lagrangian PT up to 1-loop (green) and 2-loop (red) 
results respectively, i.e., [S,N-hody{f) — CPT('^)]/'^Ar-body Note again that the error bars of the 
A^-body results are rather small, and are visually hard to see in figure 3. 

Apart from the large scales beyond the location of acoustic peak, where the reliability 
of A^-body simulations becomes subtle due to the limited simulation boxsize [6], a good 
agreement between the A^-body simulations and predictions is found for both 1-loop and 
2-loop LPT. As decreasing redshifts, the acoustic peak structure tends to be erased and the 
peak position is shifted to a small scale. These are purely the result of non-linear mode 
coupling, and the LPT explains this smearing effect quite well [37]. On the other hand, the 
shift of the baryon acoustic peak position is at a sub-percent level. Even at redshift z = [41] , 
it is difficult to quantify this effect in our results of A^-body simulations. An interesting point 
is here that the correlation function is hardly affected by the 2-loop correction compared to 
the power spectrum. This trend has been also reported in other PT treatment (e.g., closure 
theory [6] and RPT [30]). The difference between A^-body simulations and both 1-loop and 
2-loop results of LPT is small enough around the acoustic peak (r ~ W5h~^ Mpc), and the 
accuracy of both LPT predictions reaches at several percents level, which is sufficient for 




(4.3) 
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Figure 3. Two-point correlation functions given at redshifts z = 3 (top left), 2 (top right), 1 (bottom 
left) and 0.5 (bottom right). Two-point correlation functions from linear (black dotted), 1-loop LPT 
(green dashed) and 2-loop LPT (red solid) are compared with A^-body simulations (black symbol with 
error bars). 



ongoing and future surveys [6]. A closer look at the baryon acoustic peak at lower redshifts, 
however, reveals that 2-loop LPT tends to oversmear the peak structure, and the agreement 
with iV-body simulation looks somewhat degraded. Although we have not yet succeeded 
to identify the origin of this discrepancy, the result may suggests either the significance of 
higher-order corrections at lower redshift or the need for delicate numerical treatment in 
evaluating the multi-dimensional integration. By contrast, a noticeable effect of the 2-loop 
corrections can be seen on small scales (r < 50 Mpc), where the nonlinear enhancement 
of the correlation function becomes prominent at lower redshifts. Even though the 1-loop 
LPT tends to deviate from A^-body results, the 2-loop LPT reproduces the A^-body trends 
quite well within 1% precision. 

5 Summary 

In this paper, we present an improved prediction of the nonlinear perturbation theory via 
the Lagrangian picture, originally proposed by ref. [37]. Based on the previous result [37], we 
derive a general relation between the power spectrum in SPT and that in LPT for arbitrary 
loop-order. Using this relation, we explicitly write down the analytic expression of the power 
spectrum in LPT up to the 2-loop order in both real and redshift spaces. 
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Figure 4. Fractional differences between iV-body and Lagrangian PT results, [£,N-hody{r) 
CpT(r)]/Cw-body(r); 1-loop: green, 2-loop: red. 



Comparing the results in LPT with A^-body simulations in real space, precisely, we 
quantitative study the validity range of LPT at various redshifts. For power spectrum, 
the higher-order corrections enhance the power of high-fc. Including the 2-loop corrections, 
the range of validity, where the LPT prediction agrees with A^-body simulation within 1% 
precision, is improved by a factor of 1.0 {z = 0.5), 1.3 (1), 1.6 (2) and 1.8 (3), compared with 
the 1-loop LPT. On the other hand, the two-point correlation functions around the baryon 
acoustic peak are less sensitive to the higher-order corrections in all redshift range. This 
implies that the 1-loop LPT is sufficient to accurately describe the baryon acoustic peak in 
correlation functions. 

Finally, we should note remaining issues and future prospects of our work. We postponed 
the evaluation of redshift-space power spectrum and the check of its accuracy. The expression 
of redshift-space power spectrum involves a difficult term to evaluate, and we need to develop 
an efficient algorithm to compute it. On the other hand, in practice, it is inevitable to properly 
take account of the galaxy biasing in the real BAO survey. The LPT provides a way to treat 
the nonlinear galaxy biasing in a self-consistent manner, and the previous study [53] has 
reported the advantage of LPT based on the 1-loop calculations. The extension to the 2-loop 
calculations would be presumably a straightforward task, and will be reported elsewhere. 
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